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Abstract 

We study two types of generalized Baxter- Wu models, by means of transfer-matrix and Monte 
Carlo techniques. The first generalization allows for different couplings in the up- and down 
triangles, and the second generalization is to a g-state spin model with three-spin interactions. 
Both generalizations lead to self-dual models, so that the probable locations of the phase transitions 
follow. Our numerical analysis confirms that phase transitions occur at the self-dual points. For 
both generalizations of the Baxter- Wu model, the phase transitions appear to be discontinuous. 

PACS numbers: 05.50.+q, 64.60.Cn, 64.60.Fr, 75.10.Hk 
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I. INTRODUCTION 



In general, systems in the universality class of the two-dimensional 4-state Potts model 
display critical singularities that are modified by logarithmic correction factors. A satisfac- 
tory explanation of this fact is provided by the renormalization scenario due to Nienhuis 
et al. [Ij]. It explains the logarithmic factors [2! as arising from the second temperature 
field, which is marginally irrelevant. It also shows that the 4-state Potts behavior without 
logarithmic factors can only occur at special points in the parameter space, where the two 
leading temperature fields simultaneously vanish. The exactly solved Baxter- Wu model 
precisely fits such a location in parameter space: it belongs to the 4-state Potts class and its 
leading critical singularities do not have logarithmic factors. Its reduced Hamiltonian reads 

AV 

where /3 = l/lksT) is the inverse temperature, and the sum is over the up- and down 

triangles of the triangular lattice, and the site labels i, j and k refer to the three spins at 

the vertices of each triangle. Each spin assumes the Ising values ±1; this is emphasized 

by the superscript I of the coupling K^. At low temperatures, the model is in one of four 

long-range ordered phases, where most triangles have an even number of — spins. While 

the common type of interaction between spins in magnetic materials is of the two-spin type, 

three-particle interactions such as in the Baxter- Wu model have been used to describe the 

fl 

shape of face-centered cubic crystal surfaces 

This work investigates two different generalizations of the Baxter- Wu model. First we 
consider the case that the couplings in the up- and down triangles are different (see Fig. [1]), 
i.e., 

pn = -Kl ^ SiSjSk - Kl'Y SiSjSk (2) 
A V 

where the sums are over the up- and down triangles of the triangular lattice respectively. 
The introduction of another temperature-like parameter makes it likely that this model will 
have a critical line parametrized by the ratio of Kl and Kl. The fourfold degeneracy of 
the ground state persists for Kl 7^ Kl, so that it may seem plausible that the model still 
belongs to the 4-state Potts universality class. We shall attempt to provide a more definite 
judgment by means of a numerical investigation. 
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FIG. 1: A 6 X 6 triangular lattice. The lattice is divided into a honeycomb (open and filled circles) 
and a triangular sublattice (pentagons), which are dual to each other. The honeycomb lattice is 
bipartite. 

For the second generalization it is useful to write Eq. ([T]) in terms of two-state Potts 
variables cTj = (sj + 3)/2 = 1 or 2: 

pn = -Kj2^2{cT, + a, + ak) (3) 

AV 

where 52{x) = if a; is odd and 1 if x is even, and K = 2K^. The sum is over all up- 
and down triangles. Eqs. ([1]) and ([3]) differ by an additive constant that is irrelevant for 
the present purposes. It is now straightforward to generalize the model in terms of g-state 
variables with values cxj = 1, 2, ■ ■ ■ , q: 

(3n = -Kj2 + + ^-fc) (4) 
AV 

where Sq{x) = 1 if (x mod q) = and Sq{x) = otherwise. This model can also be consid- 
ered as a generalization of the g-state Potts model [Si] to 3-spin interactions, because the pair 
couplings of the original Potts model on a bipartite lattice can be written as —K6q{ai + crj)- 
But the model (jlj) does not obey the g-fold permutation symmetry Sg of the Potts model 
for general q. Its symmetry group is Zq® Zq® Z2® Sj, where the g-state clock symmetries 
Zq are generated by the operation o"j ^ cTj + 1 mod g, independently for two of the three 
sublattices, Z^ is generated by the operation cij ^ g + 1 — (Xj on all sites, and 1S3 is the 
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symmetric group of the permutations of the three sublattices. The latter symmetry results 
from the spatial symmetries of the lattice, namely reflection and translation or rotation, 
which can permute the three sublattices, while leaving (31-L invariant. 

It is obvious that the degeneracy of the ground state increases as with the number q 
of spin states, so that one may expect that the model will display a discontinuous ordering 
transition for q > 2. However, the special nature of the critical Baxter- Wu model, i.e. the 
model of Eq. (jlj) for g = 2, namely the vanishing of the marginal temperature field, opens 
the possibility of another scenario. After a mapping on the Coulomb gas Q], the marginal 
temperature field translates into the fugacity of the e = 4 electric charges. Thus the q = 2 
transition maps precisely on the point of the Gaussian fixed line where the electric charges 
are absent, and there seems to be a real possibility that this is also the case for other values 
of q. Since one expects that the Coulomb gas coupling increases with q, the electric charges, 
which are marginal at g = 2, must be relevant for q > 2, and would drive the ordering 
transition first order. But, if these charges remain absent, the transition still takes place on 
the Gaussian line, and must be critical. 

For this reason it is interesting to investigate the character of the ordering transition for 
q > 2. There are existing results due to Alcaraz et al. 0, Q] who investigated a different 
generalization of the Baxter- Wu model, namely, to a p-state clock model. For the case 
p = 3, their model is equivalent with our q = 3 model. They concluded that the transition 
is first-order for p = 3, on the basis of approximate renormalization calculations, and Monte 
Carlo calculations starting in the ordered and the disordered states, displaying changes of 
phase. 

The property of self-duality plays an important role in the present work, because knowl- 
edge of the critical point greatly facilitates the numerical analyses. Its derivation is the 
subject of Sec. Hi] where we formulate a relatively simple proof of self-duality for a class of 
models that includes both generalizations of the Baxter- Wu model mentioned above. In 
Sec. IIVI we present our numerical analysis of the q = 2 model with two different couplings, 
and in Sec. |V] we report our findings for the q = 3 and 4 models with uniform couplings. 
The conclusions of our analyses are listed in Section IVI[ 
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II. DUALITY OF g-STATE MODELS WITH MULTISPIN INTERACTIONS 



Self-duality is a useful tool to locate phase transitions. If a single phase transition occurs 
as a function of temperature, then the transition must occur at the point where the tem- 
perature variable K and the dual temperature variable K coincide. In the case of self-dual 
models with two variables Ki and K2-, the transitions tend to occur on the self-dual line in 
the Ki, K2 plane, i.e., in a point that maps onto itself under duality. 

Duality was first found for the square-lattice Ising model by Kramers and Wannier 
who correctly predicted the critical point at 

^c = ^ln(l + v^) (5) 



and since then many more derivations have been reported. Gruber et al. [10| have formu- 
lated a very general proof that includes all systems studied in the present work. For the 

convenience of the reader we shall provide a simple proof that is less general than that of 
□ 

Gruber et al. [lyj, but still more general than actually required for the models under the 
present investigation. 

Simpler, and less general versions of the proof given by Gruber et al. appear elsewhere in 
the literature. Examples are the two-dimensional Ising model with pair interactions in one 
direction and multispin interactions in the perpendicular direction (see Refs. llj, and 12 1 
for a generalization to g > 2 Potts models with similar interactions). 

The present derivation of self-duality applies to a system of g-state variables located on a 
simple hypercubic lattice. The variables are denoted dr and take the values 1, ■ ■ ■ , g. Their 
interactions are described by a Hamiltonian of the general form 

-Pn = K,Y,^, ( E ^r+a, ) + 5^ 5J ^r+b, ) (6) 

r \j=l / r \j=l J 

where r is a lattice vector and and are vectors pointing from position r to the sites of 
the variables participating in the interaction assigned to site r. There are two multiparticle 
interactions per site, one with n participating sites and another with m sites. The class 
includes the square-lattice Potts model with nearest-neighbor interactions, after a suitable 
renaming a — g — cr of the q states on one of the two sublattices. It also includes the 
Baxter- Wu model for n = m = 3, ai = —61 = (0,0), 02 = —62 = (l^O), 03 = 63 = (0, 1), 
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g = 2, and Ki = K2. The partition function for our class of models takes the form: 



{ar} r 



1 + Vi6q ^ CTr+a, 



i=l 



1 + V26q j ^ (Tr+bj 



(7) 



where = exp{Ki) — 1 and f 2 = exp(K2) — 1. Each (5g-function in Eq. ([71) can be substituted 
by its Fourier representation 



1 



±2TTizt/q 



t=l 



and each "1" in Eq. (jTl) can be replaced using the identity 

<? 

t=i 



(8) 



(9) 



The effect of these substitutions is that two new variables tr and t'j. are introduced on each 
site r, for the n- and m-particle interactions, respectively. This leads to 



V2 



exp 



2m 



(10) 



After reordering the summations and the products and collecting terms with the same a, 
we obtain 



{t,t'} 



J^(l/g)exp 



(27rzar/g) l^tr-a, -J^Cb, 



i=l 



(11) 



where N is the total number of sites in the lattice. A nice property of Eq. f[TT|) is that the 
degrees freedom on different sites r are completely independent, and thus the summation 
over the becomes very easy. Using again Fourier-transformation ([H]), one has 



tr-a.-2.t;_bj ■ (12) 

. i=l j=l 

In short, the original g- valued variable cXr has been integrated out. The price paid is the 
introduction on each site of two new g-valued variables tr,t'j. with an additional (5-function 
constraint. 



Next, one introduces a new g-state variable cXr on each site, and let t: 



tr = ^ (Jr-hj mod q, 



(13) 



which will be feasible for appropriate boundary conditions. The 6 function connecting tr 
and t'j. in Eq. (IT^ is satisfied if 



tr = ^ O-r-a, Hiod q . 



(14) 



1=1 



As the number of new variables dr is equal to the number of old variables t and t' reduced 
by the number of constraints on t and t' imposed by the rightmost S function in Eq. fll2l) . 
we expect that the dr are determined up to a trivial shift. After an inversion of the lattice, 
Eq. f|T2|) takes the form 



N 



En 



\ 1 J 



cr 



1 \ j=l 



(15) 



Comparison with Eq. ([7]) shows that Z{vi,V2) satisfies the self-duality relation 



Z{Vi,V2) 



ViV2 



N 



Z{vi,V2) 



with V1V2 = q , V2V1 = q . 



(16) 



The dual set of coupling constants {Ki,K2) obey 



vi = 6-^1-1 and V2 = e^^ -1 . 



(17) 



Each point on the line 

V1V2 = q (18) 

is mapped onto itself, and we find, for the case Vi = V2 the symmetric self-dual point as 
V = y/q or 

K = ln(l + ^) . (19) 

In this self-dual point the average number of satisfied multiparticle interactions ( "satisfied" 
means that the sum modulo q of the spins coupled by the interaction vanishes) per site, if 
unique, is found from the derivative of InZ with respect to the coupling constants at the 
self-dual point. In the case of a first-order transition on the self-dual line, this yields the 
mean of the values in the disordered phase and in the ordered phase. 
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For q = 2 models defined in terms of Ising spins Sj = ±1, one lias to take into account 
the factor 2 between the "Potts" and "Ising" couplings, as appearing under Eq. (l3])-i.e., 
V = exp{2K^) — 1. In the Ising case, the equation for the self-dual line Eq. (fT8l) may be 
written as 

smh{2Kl) smh{2Kl) = 1 (20) 

In many cases, the self-dual line, or a part of it, is the locus of a phase transition. The 
existence, uniqueness, and character of a phase transition, however, are not determined by 
self-duality. For that purpose, additional calculations are required. For several Ising models 
with multispin interactions and a field (m = 1), including three-dimensional models, Blote 
et al. [l^ found discontinuous transitions on a part of the self-dual line, with a gas-liquid 
like critical point at the end of the first-order range. For a two-dimensional system with 
pair interactions in one direction and multiparticle interactions between p particles in the 



perpendicular direction, Zhang and Yang 12| concluded, from Monte Carlo calculations. 



that a phase transition occurs at the self-dual point, and that it is first-order for all g > 2 

if n > 2. Also in the case of the n-state clock model with three-particle interactions on 

n 

the triangular lattice, Alcaraz et al. found from Monte Carlo calculations [7,] that phase 
transitions occur at the self-dual points for n = 2 and n = 3. 



III. NUMERICAL METHODS 

We investigate the generalized Baxter- Wu model ([6]) on the triangular lattice, both by 
transfer-matrix method and by Monte Carlo simulations. 

A. Transfer-matrix 

The transfer-matrix techniques used in this work are adequately described in the liter- 
ature, although the information is divided over different papers. The essential parts are 

n n n 

explained in Refs. [M], [I5j and [16[. Here we only add a few general and specific remarks 
for the convenience of the reader. From a few of the leading eigenvalues of the transfer 
matrix, one can calculate the free energies, the magnetic and energy-like correlation lengths 
of L X oo systems. For the case g = 2 we could perform such calculations up to finite sizes 
L = 27. The geometry is that of the triangular lattice wrapped on a cylinder, with one set 
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of edges perpendicular to the axis of the cyhnder. The finite size L is specified such that 
the circumference of the cyhnder is spanned by L lattice edges. 

Here we use the true triangular lattice, instead of the representation as a square lattice 
with one set of diagonal bonds, as used in Sec. [Tll Since, after adding one layer of spins, the 
lattice is shifted by a half lattice unit along the finite direction, we chose a transfer matrix 
that adds two layers of spins and applies an additional reverse shift operation, in order to 
ensure that the transfer matrix commutes with the lattice reflection as specified below. Such 
commutation relations allow one to find a common set of eigenstates of the transfer matrix 
and a symmetry operator. 

The transfer matrix acts on a vector space with vector indices representing the state of 
a row of L Ising spin variables. For q = 2, the vector indices can thus be written as binary 
numbers • • with bk = {sk + l)/2. For g = 3 one uses ternary numbers, etc., 

but here we shall use the language for binary numbers. The transfer matrix calculations 
focus on three eigenvalues, namely the largest one Aq, the "magnetic" one A^, and the 
"thermal" eigenvalue Xt- These eigenvalues are defined in the usual way, by means of the 
group of symmetry operations that leave the Hamiltonian invariant, but permute the ordered 
phases. The thermal eigenvalue, like the largest eigenvalue corresponds to an eigenvector 
fully invariant under these symmetry operations. The magnetic eigenvalue is the largest 
one with an eigenvector that changes under these symmetry operations. In this model the 
relevant symmetry group is generated by the allowed permutations of the q states, and by 
lattice symmetries that permute the three sublattices. As the transfer matrix breaks some 
of the latter symmetries, we replace the full symmetry group by the subgroup that is not 
violated by the transfer matrix. 

The analyses based on A^ and Am are similar. We proceed as follows for the case of 
Am- The magnetic correlation function gm{r) as a function of the distance r in the length 
direction of the cylinder is defined as gm{r) = (soSr). For sufficiently large r, gm{f) decays 
exponentially on a length scale ^m that depends on L and the couplings, i.e.. 




(21) 



and is determined by the eigenvalues Aq and A^ of the transfer matrix: 




(22) 
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The geometric factor y/S allows for the thickness of two layers added by the transfer matrix, 
expressed in the same unit as the finite size L. With the help of Cardy's conformal mapping 



17l | of the infinite plane on a cylinder with a circumference L, one can now, for a system at 
criticality, relate the magnetic scaling dimension Xh, which describes the algebraic decay of 
the correlation function in the infinite system, to ^m- Defining the scaled gap Xh{Ki, K2, L) 



by 



X,{K,, K,, L) ^ , (23) 

27r4^(Ai, A2,Lj 



and using finite-size scaling [l8|, one finds that, at criticahty, 



Xh{Ki, K2, L)=Xh + hiiy^ + 62^^^ + • • • (24) 

where the correction terms arise from irrelevant fields, whose presence means that 

conformal invariance applies only in the limit of large length scales. Since the irrelevant 
exponents satisfy Ui < 0, Xh{Ki, K2, L) converges to X^ with increasing L, and numerical 
estimates of Xh can be obtained from the finite-size data that can be calculated for a range 
of system sizes. 

For a system that is not critical due to the presence of some relevant scaling field, a 
term with a positive power of L appears in Eq. (!24l) . which will lead to crossover to different 
behavior, for instance described by a zero-temperature or an infinite-temperature fixed point. 
A finite-size analysis of the quantity Xh{Ki, K2, L) may thus show whether or not the system 
is critical, and if so, provide information on the universality class of the model. 

The analysis of the temperature dimension Xt from the energy-like correlation length 
similarly uses the eigenvalue A^. The calculation of this eigenvalue, with the same symmetry 



as Aq, is described in Ref. 



15| 



B. Monte Carlo algorithm 

Simulation of the generalized Baxter- Wu model on the triangular lattice can simply em- 
ploy the standard Metropolis method which involves single-spin updates only. However, a 
more efficient algorithm-a Swendsen- Wang-type cluster Monte Carlo method-can be formu- 
lated, which was already described for the Baxter- Wu model in Ref. [l9 |. 

To construct such a cluster method, one first divides the triangular lattice T into three 
sublattices Cti, and which are triangular. The union of any two sublattices form 
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a honeycomb lattice Ch which is dual to the remaining triangular lattice (see Fig. [T]). The 
partition sum of a generalized Baxter- Wu model can then be written 

{o-fe} {o-j.CTj} (ij) 

where the product is over every edge of the honeycomb sublattice Ch-, and k and k' are the 
two neighboring sites on the remaining triangular sublattice, on either side of edge {ij). The 
statistical weight associated with each edge {ij) is then 

[1 + vi5q{ai + CTj + ak)]\l + V25q{ai + aj + a^i)] = 
J2 [viSq{a, + a,+ak)]^^^ ^ h5g(cri + cr^ + a^O]^''? (26) 

where fi = exp{Ki) — 1 and f2 = exp{K2) — 1, and the convention 0° = 1 has been used. 
Thus, by introducing two bond variables bl^\blf for every edge of Ch, and replacing the 
corresponding edge weights in Eq. ( l25l) according to Eq. ( 126|) . one obtains a joint spin-bond 
model. 

The Swendsen- Wang-type cluster method can be adapted to simulate the joint spin- 
bond model. Two basic steps are involved: the bond- and the spin updates. Given a spin 
configuration, Eq. ( !26l) tells that the bond updates can be performed as in a uncorrelated 
bond percolation: the bond-occupation probability is p = fi/(l + vi) for hf^ on each edge 

/ox 

with a satisfied up triangle, and p = ^2/(1 + V2) for bij on each edge with a satisfied 
down triangle, and p = otherwise. Given a bond configuration, Eq. (!26|) tells that spin 
configurations satisfying the 6 functions have equal probability. Making use of the fact that 
the honeycomb lattice is bipartite, one can formulate the following algorithm. 
Cluster algorithm, version 1: 

1. Sublattice division. Randomly with equal probability label the three sub lattices as 1, 
2 and 3. Then merge two sublattices into a honeycomb lattice Ch = Ct2 U £t3- 

2. Bond update. On each edge (ij) of the honeycomb sublattice Ch, place an occupied 
bond with probability p = 1 — e~^^~^'^ if both the up- and the down-triangles are 
satisfied, p = 1 — e^^^ if only the up triangle is satisfied, p = 1 — e^^^ if only the down 
triangle is satisfied, and p = otherwise. 
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3. Cluster construction. A cluster is defined as a group of sites connected tlirougli occu- 
pied bonds, irrespective of colors. Decompose the lattice Ch into clusters (including 
single-site clusters). 

4. Spin update. All the spins on the triangular sublattice Ct^ are left unchanged. Ran- 
domly with uniform probability choose a value r = 0, 1, 2, ■ ■ ■ , g — 1. Independently 
for each cluster, update the spins on sublattice £^2 according to a ^ {a + t) mod q, 
and the spins on £7-3 according to a ^ {a + q — t) mod q. 

This completes one Swendsen- Wang-type cluster step, and a new spin configuration is ob- 
tained. Other choices are possible to choose r in step 4, for instance r = with probability 
1/2 and the other values of r with probability 1 / (2g — 2). The choice r = with probability 
and the other values with probability 1/g is only applicable for q > 2. 

For the special case q = 2 the cluster algorithm can be made more efficient. Conditional 
on the frozen spin configuration on sublattice Ct^, the honeycomb sublattice of the q = 2 
generalized Baxter- Wu model reduces to an Ising model with position-dependent couplings 
on the honeycomb lattice Ch- 

pn\{^^}^keCT^ = - X] SiSj{Klsk + Klsk') = - KijSiSj {s = ±1) (27) 

where the meaning of k and k' is the same as in Eq. fl25|) . The effective coupling Kij is defined 
by the right-hand side of this equation, and can be ferromagnetic or antiferromagnetic, 
depending on the spin variables Sk and Sf. On the basis of Eq. ( 127|) . the ''bond-update'^ step 
can be reformulated as follows. 
Cluster algorithm, version 2: 

2. Bond-update. On each edge {ij) of Ch, place an occupied bond with probability 
p = max[0, 1 — exp{—2siSjKij)]. 

The other steps are equal to those of version 1. An occupied bond can be either "ferromag- 
netic" or "antiferromagnetic" (between spins of opposite signs). A cluster in version 1 may 
be further decomposed into several clusters in version 2. 

We found that version 2 performs much better than the Metropolis algorithm, in the 
sense that a simulation using the cluster method yields statistically more accurate results in 
a given time. For the q = 2 case with Ki = K2, we found the dynamic exponent z as about 
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1.1, which is close to the Li-Sokal bound 



20( 1 z > 2yt — 1 = 1. For the self-dual points with 



Ki 7^ as well as those with g > 2, a further increase of the slowing down was observed. 

We mention that a single-cluster version of the algorithm can also be formulated. How- 
ever, we found that it does not further improve the efficiency. In fact, for the q = 2 case with 
Ki = K2, the dynamic exponent appears to exceed that of the full cluster-decomposition 
method. 



IV. RESULTS FOR q = 2 AND Ki / K2 

For the present case g = 2 we use the Ising notation for the condition of self-duality 
as expressed by Eq. (l20l) . Our numerical analysis divides into two parts. The transfer- 
matrix results are described in subsection IIV A[ The Monte Carlo investigation is reported 
in subsection IIVBI 



A. Transfer-matrix results 



We calculated the scaled gaps at the self-dual points with Kl = K\, and K\ = 0.5, 0.6, 
■ ■ ■ , 1.2, for system sizes up to L = 27. The system sizes were restricted to multiples of 
3, because otherwise three of the four ground states do not fit in a lattice with period L. 
For the pure Baxter- Wu model at criticality, with Kl = Kl = [ln(l -|- v^)]/2, we find that 
the finite-size data for the scaled gaps rapidly approach the exact values X^, = 1/8 and 
Xt = 1/2. Three-point fits according to 

Xf,{L)c:^Xf, + aLP (28) 



followed by iterated fits as described in Ref. [15|] reproduce the exact values up to about 



10^^ For Kl Kl, the finite -size dependence of the scaled gaps becomes stronger while 
the signs of convergence disappear, at least for a certain range of Kl/ Kl. This is illustrated 
by the finite-size data in Table [H 

These data show that the scaled gaps for the larger system sizes tend to move away 
from the exact values for the Baxter- Wu model when Kl increases. Moreover, the finite-size 
dependence, as indicated by the difference of the scaled gaps for the two largest system sizes, 
increases with Kl, except for the entry for Xt at largest value Kl = 1.2. Another significant 
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L 

FIG. 2: (Color online) Scaled gaps a function of system size L, for various pairs of couplings 
{K\,K2) proportional to the self-dual pair with K\ = I (and K\ = 0.13617- ••)• From top to 
bottom, the data apply to 0.95, 0.975, 1, 1.025, and 1.05 times the self-dual couplings. These data 
suggest the presence of a phase transition at or near the self-dual point. 

phenomenon is that the exponent p obtained by the three-point fit for the largest available 
system is positive for a range of K\, i.e., there are no longer signs of convergence with L. 
Only in the case of i^'} = 1.2 the exponent becomes negative at the largest available system 
size, which is a sign that the renormalized system is approaching an attractive fixed point. 
The presence of a phase transition can be deduced from the scaling behavior of the scaled 
function of temperature. The scaled magnetic gaps were calculated at couplings 
equal to 0.95, 0.975, 1, 1.025, and 1.05 times the self-dual pair {K\, Kl) with K\ = 1.0. These 
data are shown in Fig. [2J For the smallest coupling and largest values of L the behavior 
tends to become linear as a function of L, which corresponds with a correlation length 
that becomes constant, as expected in a disordered phase. For the largest couplings, the 
scaled gap tends rapidly to zero, which corresponds with a long-range ordered phase. This 
crossover with increasing L, which is to the high temperature phase or to the ordered phase 
for {K\, KI) smaller or larger than the self-dual pair respectively, confirms the presence of a 
phase transition at the self-dual coupling. 



0.950 □ 
0.975 



1.000 A 
1.025 V 
1.050 
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TABLE I: Results of transfer-matrix calculations of the scaled magnetic (left hand side) and energy- 
like (right hand side) gaps at the self-dual points of the generalized Baxter- Wu model for several 
values of KI. The second column indicates the type s of the scaled gap, h for and t for Xf. 
The third column shows the scaled gap for the largest available system size L = 27. Its finite-size 
dependence is indicated in the fourth column as the difference between the scaled gaps for the two 
largest finite sizes. The effective exponent p describing the finite-size dependence of the scaled gap 
is listed in the rightmost column, based on the scaled gaps for L = 21, 24, and 27. Positive values 
of p mean that the system is renormalizing away from a fixed point. The values of X^ and Xt in 
the first line in this table are close to the exact values of the scaling dimensions; the other entries 
for Xfi and Xt have no physical meaning except describing the crossover to another fixed point, 
possibly with X^ = Xt = 0. 
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-0.23 


0.8 


0.103767 


-0.0013 


0.45 


0.324382 


-0.0044 


0.10 


0.9 


0.088070 


-0.0017 


0.62 


0.244713 


-0.0080 


0.24 


1.0 


0.068388 


-0.0033 


0.83 


0.170464 


-0.011 


0.21 


1.1 


0.048182 


-0.0046 


0.42 


0.110422 


-0.013 


0.00 


1.2 


0.031335 


-0.0051 


0.00 


0.067915 


-0.012 


-0.37 



B. Monte Carlo results 

The evidence that the symmetric Baxter- Wu model {Ki = K2, q = 2) undergoes a 
second-order phase transition is very solid from the exact solution, an exact mapping to the 



0(2) loop model on the honeycomb lattice 2l|], and the existing numerical data. 



Using the aforementioned Swendsen- Wang- type cluster algorithm (version 2), we simu- 
lated the (? = 2 generalized Baxter- Wu model at the self-dual line with Kl = 0.6 and 0.8. 
The linear system size L was taken as multiples of 6 in the range 6 < L < 192; periodic 
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boundary conditions were imposed. Several quantities were sampled, including the number 
of satisfied up (down) triangles per site —E^ {—E^), the energy density E, the specific heat 
C = L'^{{E'^) — (E)"^), and the squared magnetization, defined in analogy with the rip-state 
Potts model as 

Tip — 1 Tip 



-1 j_ 

3lE (29) 



rrip = — 

1=1 j=t+i 



where we have divided the satisfied triangles into np = groups according to the associated 
ground states, and pi, with i = l,np, is the density of triangles in the ith ground state. 
We fitted the C data by 

C{L) = a + bL^-^""^ (30) 

and the m| data by 

m|(L) = L-2^'^(a + 6L-i) , (31) 

where a and b are unknown constants. The fits yield Xt = 0.43 (2) and Xh = 0.1208 (6) for 
Kl = 0.6, and Xt = 0.30 (3) and Xh = 0.110 (2) for Kl = 0.8. The results are compatible 
with those in Table [H 

The probability distributions P for the sampled quantities are also analyzed. The distri- 
bution -P(-E'u) of the density —E^^ of the satisfied up-triangles appears to be clearly bimodal, 
but the two peaks have unequal heights. The reweighted distributions Pr were obtained 
by multiplication of -P(-Eu) with a factor e'*+''^", with a and b chosen such that Pr(-E'u) is 
normalized to 1 and that its two peaks have equal heights. This transformation takes away 
an overall gradient in the energy distribution so that the signature of a first order transition 
is clearly visible. Figure [3] shows as a function of E^, and the distance AE^ between its 
two maxima. 

For first-order transitions, we expect the following behavior of the reweighted energy 
distribution: 

and the local 



1. The difference between the maximum probability density max[P(£'u 
minimum min[P(£'u)] between both maxima increases as L increases 

2. The distance AE^ approaches to a nonzero value when L ^ oo. 

The data shown in Fig. [3] are in agreement with these conditions. The horizontal scale is 
chosen as L~^/^ because AE^ then behaves approximately linearly in the pertinent range 
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FIG. 3: (Color online) Reweighted probability distribution Pj-^E^), and the distance between the 
two-peak positions AE'u for K\ = 0.8. The height of the peaks increases with system size. 

24 < L < 384. For larger L we expect a faster type of convergence, which means that the 
extrapolation in Fig. [3] may slightly underestimate the energy discontinuity for L —>■ oo. We 
also sampled the probability distribution of the magnetization-like quantity rrip, and found 
the same type of behavior, in agreement with both conditions. In short, the evidence shown 
in Fig. [3] for the generalized q = 2 Baxter- Wu model with Kl = 0.8 is just as expected for a 
first-order transition. 
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FIG. 4: (Color online) Hysteresis loop for the squared magnetization rrip of an L = 576 system 
with K1/K2 = 5. The horizontal scale shows the couplings in units of the self-dual couplings for 
this ratio. Each data point represents a simulation of 5 x 10^ Metropolis sweeps. The results for 
increasing couplings are shown as A, for decreasing couplings as V. The lines are added for visual 
aid only. 

In the case of a first-order transition, we also expect metastable phases in a temperature 
range about the self-dual point, with lifetimes that are much larger than the time scale 
describing the jump from a metastable to a stable branch. We checked for such hysteresis 
in the model with K1/K2 = 5 by simulations sweeping slowly over ranges of couplings 
including the self-dual point. To find clear hysteresis loops, one has to simulate rather large 
systems. Results for L = 576, with data points representing simulations of a half million 
Metropolis sweeps, separated by steps of 10"'' times the self-dual coupling, are shown in 
Fig. m The hysteresis loop covers only 10~^ of the K/K^d scale, where Ksd denotes the 
self-dual couplings. 
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TABLE II: Results of transfer-matrix calculations of the scaled magnetic gaps for the q = 3 and 
4 generalized Baxter- Wu models, at the respective self-dual points with Ki = K2. The columns 
under "p" show the exponent obtained from the three-point fits described in the text. 





g = 3 


q = 


4 


L 


Xh{L) 


P 


Xh{L) 


P 


3 


0.129163 




0.13050 




6 


0.117738 


1.19 


0.10381 


1.04 


9 


0.105105 


0.71 


0.07655 


0.37 


12 


0.093650 


0.62 


0.05460 




15 


0.083255 


0.54 






18 


0.073778 









V. RESULTS FOR q>2 

A. Transfer-matrix calculations 

We have constructed transfer-matrix algorithms for the g = 3 and 4 generalization of the 
Baxter- Wu model with Ki = K2. The program is rather similar to that for the Baxter- 
Wu model, the main difference is that we have to use ternary or quaternary numbers to 
characterize a row of site variables, instead of binary numbers. As a consequence, a smaller 
range of system sizes can be handled. The finite-size data are here restricted to L < 18 for 
g = 3 and L < 12 for g = 4. 

We computed the largest eigenvalue of the transfer matrix, as well as the magnetic eigen- 
value, characterized by the antisymmetry under a lattice reflection of the corresponding 
eigenstate. Next, the correlation length and the scaled gap were obtained from Eqs. (122|) 
and ( l23l) . The results for the scaled gap are shown in Table HTl 

The behavior of the scaled gaps does not suggest convergence with increasing L. Three- 
point fits according to Eq. ( 128|) yield positive values of the exponent p. This does not agree 
well with the description of the finite-size data in terms of an attractive critical fixed point. 
It rather suggests crossover to some other, sufficiently remote fixed point. That may well be 
a discontinuity fixed point [23]. Both for g = 3 and 4, the behavior of the scaled gaps as a 
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FIG. 5: (Color online) Scaled gaps of the q = 3 generalized Baxter- Wu model as a function of 
system size L, for five different couplings in the vicinity of the self-dual point. From top to bottom, 
the data apply to 0.9, 0.95, 1, 1.05, and 1.1 times the self-dual coupling. These data suggest that 
a phase transition takes place near the self-dual point. 

function of L is similar to that found in Sec. lIVAI at intermediate values of Kj. 

Transfer-matrix calculations at couplings with Ki = K2 in the vicinity of the self-dual 
value show clear signs of transitions. The scaled magnetic gaps shown in Figs. [5] and [6] for 
q = 3 and 4 respectively, display the same type of transition behavior as found ins Sec. IIV Al 
for a g = 2 model: for couplings exceeding the self-dual value the scaled gaps tend to zero, 
and at the high-temperature side the scaled gaps are increasing with system size. 



B. Monte Carlo results 



Also in this case we employ Monte Carlo simulations to obtain independent and additional 
evidence about the character of the phase transitions. In addition to the evidence already 
reported by Alcaraz et al. 7|, Isl], it remains to be investigated whether hysteresis is present, 
and whether one can extrapolate the energy discontinuity to the thermodynamic limit. 

We employed the Metropolis method as well as the cluster algorithm defined in Sec. IIII B[ 
However, in the present case q > 2, the efficiency of the cluster method is not much different 
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FIG. 6: (Color online) Scaled gaps of the q = A generalized Baxter- Wu model as a function of 
system size L, for five different couplings in the vicinity of the self-dual point. From top to bottom, 
the data apply to 0.9, 0.95, 1, 1.05, and 1.1 times the self-dual coupling. These results are similar 
to those for the q = 2> model, but the gap at the self-dual point decays more rapidly with L for 
g = 4. 

from that of the Metropolis algorithm. 

We first simulated the Ki = K2 self-dual point of the g = 3 model, and sampled the 
energy distribution for a number of system sizes that are multiples of 3. The energy E is 
defined as minus the density of satisfied triangles per site. Again the distribution has two 
unequal peaks, but their separation is wider than in the q = 2 case. The reweighting was 
done by multiplication of the histogram with e"'*"''^. The reweighted distribution Pr{E) is 
shown in Fig. [7] for several system sizes. The local minimum between the peaks decreases as 
a function of L. In the range of finite sizes covered by our simulations, the distance between 
the peaks approaches a nonzero constant approximately as 1/L, as shown in Fig. [HI Such 



behavior was also found by Lee and Kosterlitz [22| for the first-order transition of the g > 4 
Potts model. The average of the two peaks, also shown in this figure, extrapolates within 
numerical uncertainty to the value 1 + 1/ predicted by self-duality. 

Next, we performed similar simulations of the g = 4 model at the self-dual point. The 
reweighted probability distribution Pr{E) is shown in in Fig. [H] for several system sizes. The 
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FIG. 7: (Color online) Reweighted probability distribution Pr{E) for the (7 = 3 model. Data are 
shown for system sizes L = 6, 12, 24 and 48. Data points for the same system size are connected 
by a curve for the purpose of clarity. The heights of the peaks increase with system size. 
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I ' ' ' 1 
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FIG. 8: (Color online) Distance along the energy scale between the peaks (A), and minus the mean 
of the peaks (□) of the energy histogram of the g = 3 model versus inverse system size. Duality 
predicts the mean of the peaks at the position marked by O. 
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FIG. 9: (Color online) Reweighted probability distribution Pr{E) for the q = A model. Data are 
shown for system sizes L = 6, 12 and 24. Data points for the same system size are connected by a 
curve for the purpose of clarity. The height of the peaks increases with system size. 

distances between the maxima of the histogram are shown in Fig. [10] as a function of the 
inverse system size. They extrapolate to a nonzero constant. The average peak positions, 
also shown in Fig. [HI agree well with the value 3/2 predicted by duality. Also these data 
agree with the expectations for a first-order transition, and even more strongly so than in the 
g = 3 case, for instance, because the distances between the peaks of the energy histograms 
are larger. 

To test for the presence of hysteresis, we performed Monte Carlo simulations of the g = 3 
and 4 models, varying the temperature in a region close to the symmetric self-dual point. 
Each data point involved a simulation of 2 x 10^ Metropolis sweeps, of which the first 10^ 
were used for equilibration. The results for the magnetization-type quantity m| are shown 
in Figs. [Uland Fig. [121 They display a small hysteresis loop for g = 3, covering only a half 
percent of the K scale, and stronger hysteresis effects for g = 4. 
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FIG. 10: (Color online) Distance along the energy scale between the peaks (A), and minus the 
mean of the peaks (□) of the energy histogram of the q = 4 model versus inverse system size. 
Duality predicts the mean of the peaks at the position marked by O. 

VI. CONCLUSION 



The numerical results presented in Sec. IIV Al for the Baxter- Wu model (g = 2, Ki = K2) 
clearly converge to the known exact values Xt = 1/2 and Xh = 1/8. For Ki 7^ K2 deviations 
from this behavior are observed, and the dependence of these estimates on the finite size 
L is considerable when Ki and K2 are sufficiently different. At first sight, this situation 
may seem similar to the poor convergence observed for some models in the 4-state Potts 
universality class, see e.g. Ref. [l^. 

However, there are also significant differences. First we note that, except for ratios K1/K2 
close to 1, the differences in the finite-size estimates for Xt and Xh tend to increase with 
increasing system size. Second, the finite-size estimates for Xt and Xh are smaller than the 



exact values 
model 
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25|. 



'or the Baxter- Wu model, instead of larger as observed for the g = 4 Potts 



The interpretation of these observations is suggested by the renormalization flow diagram 
for the surface of phase transitions of the dilute two-dimensional Potts model proposed by 
Nienhuis et al. [l]. The parameter space of that work involved the chemical potential v of 
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FIG. 11: (Color online) Hysteresis loop of the magnetization- like quantity rn-p for the g = 3 model 
with size 60^. The horizontal scale shows the coupling in units of the self-dual coupling. The 
results for increasing couplings are shown as A, for decreasing couplings as V. The lines are added 
for visual aid only. 

vacant sites and the number of Potts states q. The mapping of the Potts model onto the 



random-cluster model 26|] enables one to treat g as a continuous variable. Since vacant sites 
in the Potts model are dual to multisite interactions [27], the parameter v may as well be 
interpreted as a scaling field depending on the type of interactions. At g = 4, the field v 
becomes marginal ^ at the critical point. 

We reproduce this flow diagram adapted to our purposes, in Fig. [131 The g = 4 
Potts model is located at a value of v smaller than that at the g = 4 fixed point, and is still 
attracted by it, although marginally. This explains the slow finite-size convergence, and the 
logarithmic factors of the g = 4 Potts model. The Baxter- Wu model is located at the g = 4 
fixed point. 

The introduction of a difference between K\ and such that the condition of self- 
duality is still satisfied, allows for the possibility that the location of the model in Fig. [T3l 
changes. The coordinate g will remain unchanged, but a priori there does not seem to be a 
way to tell whether the model will move up or down in the diagram, or perhaps will keep 
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FIG. 12: Hysteresis loop of the magnetization-like quantity nip for the g = 4 model with size 
60^. The horizontal scale shows the couplings in units of the self-dual coupling. The results for 
increasing couplings are shown as A, for decreasing couplings as V. The lines are added for visual 
aid only. 

its location. But, since the finite-size estimates of and Xt for the Kl ^ K\ models and 
those for the g = 4 Potts model lie on opposite sides with respect to the Baxter- Wu model, 
we may locate the K\ ^ K\ models at a value of v exceeding that of the Baxter- Wu model. 



as indicated by "2C" in Fig. [131 Therefore they flow to the discontinuity fixed point 23 1 
located at large f , so that the phase transition is discontinuous. In view of the symmetry 
between Ki and K2., the marginally relevant field v can, in lowest order, not depend linearly 
on Kl — K2 near the 4-state Potts fixed point, and one expects a contribution as {K\ — K}^'^. 
This is consistent with the very weak dependence of the finite-size data in Table [T] on small 
differences K\ — K\. 

Thus we conclude that the generalized Baxter- Wu model with different couplings de- 
scribed by the Hamiltonian ([2]) undergoes a phase transition at the self-dual line for q > 2, 
and that the phase transition is first order for Ki 7^ K2, although extremely weakly so when 
the difference Ki — K2 is small. Even for a rather large difference Ki/ K2 = 5, we find (see 
Fig. HI) a very narrow hysteresis loop. 
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FIG. 13: (Color online) Renormalization flow in the plane of phase transitions of the dilute random- 
cluster model, parametrized by the number of Potts states q and the fugacity v of vacancies, ac- 
cording to Ref. \x\. The curve represents a line of fixed points. Its lower branch is attractive 
and describes the critical random-cluster model. The upper branch of fixed points is repulsive 
and describes the tricritical random-cluster model. When v exceeds its tricritical value, the renor- 
malization flow leads to a discontinuity fixed point, corresponding with a first-order transition. 
The position of the Baxter- Wu model (BW), of the q = A Potts model (4P) and the presently 
investigated self-dual q = 2 models with Ki ^ K2 (2C) are sketched. 



Furthermore, for q = 3 and 4 the transition is also discontinuous. This result disproves 
the possibility mentioned in Sec. [T] that the q > 2 self-dual generalized Baxter- Wu models 
renormalize to a Coulomb gas in which the fugacity of the electric charges vanishes, in which 
case algebraic critical behavior would occur. Apparently the fugacity is nonzero, and, since 
the electric charges are relevant for q > 2, the models renormalize away from the Gaussian 
line to a discontinuity fixed point. 

The first-order character of the g = 4 model, as expressed, for instance, by the energy 
discontinuity, is stronger than that of the g = 3 model. We expect the first-order character 
to grow even stronger with a further increase of q and/or the introduction of an asymmetry 
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